Homework 3

    Describe the setup and each step in your solutions with words and clearly label your final answers. Use Matlab for plotting and programming and include your code as an appendix to your problem set.

Table of Contents

▶︎
all
running...

Problem 1

    Leapfrog for scalar advection.

    (a) Repeat the numerical solution of the scalar advection equation from Problem Set #2 using leapfrog for time advancement and central differencing in space. Use the same pulse of dye and other problem parameters (just modify your old Matlab code). Plot the results for Cr=0.5\text{Cr} = 0.5 on the same c,xc, x axis at t=0,1,2 st = 0, 1, 2 \, \text{s} and describe your plot. Note: you will have to use a different method to advance the first time step. Make sure the scheme you use is stable.

    (b) Use modified equation analysis to examine the accuracy of the leapfrog time stepping scheme used above. Make sure to express the error terms using spatial derivatives. Describe the influence of the leading error terms (are they dissipative? dispersive?) on the behavior of the numerical solution compared to the exact solution.

    (c) Use Von Neumann stability analysis to determine the stability properties of this scheme. You will get a quadratic expression for σ\sigma which you can solve using the quadratic formula.

    (d) To control the oscillations in the leapfrog-central scheme, various approaches can be used. One is to use what’s commonly known as the Asselin filter (also the Robert filter) to replace cinc^n_i with

cin+α2(cin+1−2cin+cin−1)c_i^n + \frac{\alpha}{2}\left( c_{i}^{n+1} - 2c_i^{n} + c_{i}^{n-1} \right)

at each time step. Try this for your code from problem 1 above using Asselin coefficients α=0.05\alpha = 0.05 (recommended value) and α=1.0\alpha = 1.0. What happens to your solution? Plot your results.

    本节要求解一维平流问题 (Eq. (2.23) in Hw2)

{ct+ucx=0,t≥0,  x∈[0,L],(I.C.)c∣t=0=c0(x),(B.C.)c∣x=0=c∣x=L=0,\begin{equation*} \tag{1.1} \left\{ \begin{aligned} & c_t + u c_x = 0, \quad t \ge 0, \; x \in [0, L], \\ & \text{(I.C.)} \quad \left. c \right|_{t = 0} = c^0(x),\\ & \text{(B.C.)} \quad \left. c \right|_{x = 0} = \left. c \right|_{x = L} = 0, \end{aligned} \right. \end{equation*}

式中常数 u>0u > 0. 在上一章, 此问题用前向 Euler 迎风格式 (Eq. (2.26) in Hw2) 离散. 在这里, 我们试试 Leap-frog central difference 格式,

{cjn+1=cjn−1−Cr (cj+1n−cj−1n),n∈N∗,  j∈Z∩[0,L/h],(I.C.)cj0={1,jh∈(0,0.5],0,jh∉(0,0.5],(B.C.)cjn=0,jh∈(−∞,0]∪[L,+∞),\begin{equation*} \tag{1.2} \left\{ \begin{aligned} & c^{n+1}_j = c^{n-1}_{j} - \text{Cr} \, (c_{j+1}^{n} - c_{j-1}^{n}), \quad n \in \mathbb{N^*}, \; j \in \mathbb{Z} \cap [0, L / h], \\ & (\text{I.C.}) \quad c_j^0 = \left\{ \begin{aligned} & 1, \quad & jh \in (0, 0.5], \\ & 0, \quad & jh \notin (0, 0.5], \end{aligned} \right. \\ & (\text{B.C.}) \quad c_j^n = 0, \quad jh \in (-\infty, 0] \cup [L, +\infty), \end{aligned} \right. \end{equation*}

其中 Courant number

Cr:=uτ/h.\begin{equation*} \text{Cr} := u \tau / h. \end{equation*}

    (a) 用 Matlab 实现格式 (1.2), 首个时间步的积分方案采用显式 Euler 迎风格式. 结果表明, 当 Cr=0.5\mathrm{Cr} = 0.5, 格式 (1.2) 稳定, 但出现虚假的振荡; 当 Cr=1\mathrm{Cr} = 1, 格式 (1.2) 能还原精确解, 但 Cr\mathrm{Cr} 只要略增大, 格式就失去稳定性; 只要略减小, 数值解就立即出现明显的虚假振荡.

Leap-frog central diff. 格式积分一维平流问题

    (b) 满足 Leap-frog central diff. 格式 (Eq. (1.2))

cjn+1−cjn−12τ+ucj+1n−cj−1n2h=0\begin{equation*} \tag{1.3} \frac{c_{j}^{n+1} - c_{j}^{n-1}}{2 \tau} + u \frac{c^n_{j+1} - c^n_{j-1}}{2 h} = 0 \end{equation*}

的 c~∈C6(R)\tilde{c} \in C^6(\mathbb{R}) 成立

c~t+uc~x=−τ26c~t3−h26uc~x3+O(τ4+h4),\begin{equation*} \tag{1.4a} \tilde{c}_t + u \tilde{c}_x = -\frac{\tau^2}{6} \tilde{c}_{t^3} - \frac{h^2}{6} u \tilde{c}_{x^3} + \mathcal{O}(\tau^4 + h^4), \end{equation*}

这就是格式 (1.2) 的 modified equation. 与精确方程 (1.1) 比较, (1.4a) 右端有截断误差的意义, 所以说格式 (1.2) 具有二阶时间精度和二阶空间精度. 想把 (1.4a) 右端的 c~t3\tilde{c}_{t^3} 项写成空间导数, 以明晰物理意义. 为此, 让 (1.4a) 分别被算子 ∂t2,∂tx,∂x2\partial_{t^2}, \partial_{tx}, \partial_{x^2} 作用, 得

c~t3+uc~t2x=O(τ2+h2),c~t2x+uc~tx2=O(τ2+h2),c~tx2+uc~x3=O(τ2+h2).\begin{align*} \tag{1.5} \tilde{c}_{t^3} + u \tilde{c}_{t^2 x} & = \mathcal{O}(\tau^2 + h^2), \\ \tag{1.6} \tilde{c}_{t^2 x} + u \tilde{c}_{tx^2} & = \mathcal{O}(\tau^2 + h^2), \\ \tag{1.7} \tilde{c}_{t x^2} + u \tilde{c}_{x^3} & = \mathcal{O}(\tau^2 + h^2). \end{align*}

三式联立得

c~t3=−u3c~x3+O(τ2+h2),\begin{equation*} \tag{1.8} \tilde{c}_{t^3} = -u^3 \tilde{c}_{x^3} + \mathcal{O}(\tau^2 + h^2), \end{equation*}

进而把 (1.4a) 改写为

c~t+uc~x=h26u(Cr2−1)c~x3+O(τ4+τ2h2+h4).\begin{equation*} \tag{1.4b} \tilde{c}_t + u \tilde{c}_{x} = \frac{h^2}{6} u (\text{Cr}^2 - 1) \tilde{c}_{x^3} + \mathcal{O}(\tau^4 + \tau^2 h^2 + h^4). \end{equation*}

    为获得关于 modified equation (1.4b) 的截断误差主部 (即右端第一项) 对数值解的影响的直观, 考虑波数为 ξ\xi 的单色平面波

f(x,t)=f^(t)exp⁡(iξx)\begin{equation*} \tag{1.9} f(x, t) = \hat{f}(t) \exp{(\mathrm{i} \xi x)} \end{equation*}

在系统

ft=afx3\begin{equation*} \tag{1.10} f_t = a f_{x^3} \end{equation*}

中的行为. 由满足 (1.10) 的 (1.9) 所具有的形式

f=f^(0)exp⁡[iξ(x−aξ2t)]\begin{equation*} \tag{1.11} f = \hat{f}(0) \exp{\left[ \mathrm{i} \xi \left(x - a \xi^2 t \right) \right]} \end{equation*}

看出, Leap-frog central diff. 格式 (1.3) 的误差主部 (来自 (1.4b))

h26u(Cr2−1)c~x3\frac{h^2}{6} u (\text{Cr}^2 - 1) \tilde{c}_{x^3}

会引入数值色散 (numerical dispersion), 即相速度依赖波数, 且高频波传播更快, 使初值在积分过程中变形, 可能出现高频“噪声”. 而纯平流问题 (1.1) 的精确行为, 是初值保持原样移向下游, 直至在区域边界被完全吸收.

    (c) 引入

Cjn:=[cjncjn−1]T,\begin{equation*} \tag{1.12} \boldsymbol{C}_j^{n} := \begin{bmatrix} c_j^n & c_{j}^{n-1} \end{bmatrix}^{\text{T}}, \end{equation*}

则 Leap frog central difference 格式 (1.3) 写成

Cjn+1=[−Cr0]Cj+1n+[11]Cjn+[Cr0]Cj−1n.\begin{equation*} \tag{1.13} \boldsymbol{C}_j^{n+1} = \begin{bmatrix} -\text{Cr} & \\ & 0\\ \end{bmatrix} \boldsymbol{C}_{j+1}^{n} + \begin{bmatrix} & 1 \\ 1 & \end{bmatrix} \boldsymbol{C}_j^n + \begin{bmatrix} \text{Cr} & \\ & 0 \end{bmatrix} \boldsymbol{C}_{j-1}^n. \end{equation*}

为讨论格式 (1.3) 或 (1.13) 的稳定性, 考虑波数为 ξ\xi 的单色平面波

Cjn=C^(tn)exp⁡(iξxj)\begin{equation*} \tag{1.14} \boldsymbol{C}_j^n = \hat{\boldsymbol{C}}(t_n) \exp{(\mathrm{i}\xi x_j)} \end{equation*}

在 (1.13) 积分过程中的表现, 发现有

Cjn+1=GCjn,\begin{equation*} \tag{1.15} \boldsymbol{C}_j^{n+1} = \boldsymbol{G} \boldsymbol{C}_j^{n}, \end{equation*}

其中复增长矩阵

G=G(ξ,τ,h):=[−i 2 Crsin⁡(ξh)110]\begin{equation*} \tag{1.16} \boldsymbol{G} = \boldsymbol{G}(\xi, \tau, h) := \begin{bmatrix} -\mathrm{i}\, 2 \, \text{Cr} \sin{(\xi h)} & 1 \\ 1 & 0 \end{bmatrix} \end{equation*}

的特征值 λ\lambda 满足

{λ1+λ2=trG=−i 2 Crsin⁡(ξh),λ1λ2=det⁡G=−1,\begin{equation*} \tag{1.17a} \left\{ \begin{aligned} & \lambda_1 + \lambda_2 = \mathop{\text{tr}}{\boldsymbol{G}} = - \mathrm{i}\, 2 \, \text{Cr} \sin{(\xi h)}, \\ & \lambda_1 \lambda_2 = \det{\boldsymbol{G}} = -1, \end{aligned} \right. \end{equation*}

即是特征方程

0=det⁡(λI−G)=λ2+i 2λ Crsin⁡(ξh)−1\begin{equation*} \tag{1.17b} 0 = \det{(\lambda \boldsymbol{I} - \boldsymbol{G})} = \lambda^2 + \mathrm{i} \, 2 \lambda \, \text{Cr} \sin{(\xi h)} - 1 \end{equation*}

在复数域上的两个根.

    常系数线性发展方程格式稳定的充要条件 (Eq. (3.20) in Durran[1] (2012)) 是 ∃τ0, K, ∀τ≤τ0, nτ≤T, ∣ξh∣≤π:\exist \tau_0, \, K, \, \forall \tau \le \tau_0, \, n \tau \le T, \, |\xi h| \le \pi:

∥G(ξ,τ,h)∥n≤K.\begin{equation*} \tag{1.18} \| \boldsymbol{G}(\xi, \tau, h) \|^n \le K. \end{equation*}

Eq. (1.18) 蕴含 (ρ(G)\rho(G) 是 G\boldsymbol{G} 的谱半径)

(ρ(G))n=ρ(Gn)≤∥G∥n≤K,(∀nτ≤T),⇒ ρ(G)≤K1/N=Kτ/T≤1+Mτ,(∃M=M(τ0,K), ∀τ≤τ0).\begin{align*} & (\rho(G))^n = \rho(G^n) \le \|G\|^n \le K, \quad (\forall n \tau \le T), \\ \Rightarrow & \, \rho(G) \le K^{1/N} = K^{\tau/T} \le 1 + M \tau, \quad (\exist M = M(\tau_0, K), \, \forall \tau \le \tau_0). \tag{1.19} \end{align*}

Eq. (1.19) 被称为 von Neumann 条件, 是常系数线性发展方程格式稳定的一个必要条件. 此外, 若复增长矩阵 GG 满足

∥G(ξ,τ,h)∥≤1+Mτ,(∃M>0, ∀∣ξh∣≤π),\begin{equation*} \tag{1.20} \| G(\xi, \tau, h) \| \le 1 + M \tau, \quad (\exist M > 0, \, \forall |\xi h| \le \pi), \end{equation*}

则 (1.18) holds, 故 (1.20) 是格式稳定的一个充分条件. 故若 ρ(G)\rho(G) 是 GG 的一种矩阵范数 (例如, 当 GG 是正规矩阵, 或是标量), 则 von Neumann 条件 (1.19) 成为格式稳定的充要条件. 简单起见, 且考虑物理意义 (波动不增长), 以后我们将稳定性 (充要) 条件取为 Eq. (3.21) in Durran[1:1] (2012),

∥G∥≤1,\begin{equation*} \tag{1.21} \| G \| \le 1, \end{equation*}

使 von Neumann 条件 (1.19) 加强为, 复增长矩阵 GG 的所有特征值 λ\lambda 均满足

∣λ∣≤1.\begin{equation*} \tag{1.22} |\lambda| \le 1. \end{equation*}

    现在证明, 对于 Leap frog central difference 格式 (1.13), 稳定性条件 (1.22) 等价于

∣u∣τ/h=∣Cr∣≤1.\begin{equation*} \tag{1.23} |u| \tau / h = |\text{Cr}| \le 1. \end{equation*}

先证必要性. 若 (1.22) 成立, 则由 (1.17a) 第二式知 λ1,λ2\lambda_1, \lambda_2 落在复平面的单位圆上 (即 ∣λ1∣=∣λ2∣=1|\lambda_1| = |\lambda_2| = 1); 再由 (1.17a) 第一式得

2 ∣Cr sin⁡(ξh)∣=∣λ1+λ2∣≤∣λ1∣+∣λ2∣=2\begin{equation*} \tag{1.24} 2 \, |\text{Cr} \, \sin{(\xi h)}| = | \lambda_1 + \lambda_2| \le |\lambda_1| + |\lambda_2| = 2 \end{equation*}

对一切 ∣ξh∣≤π|\xi h| \le \pi 成立, 故有 (1.23). 必要性证毕. 再证充分性. 若 (1.23) 成立, 则从 (1.17a) 看出, 对任意 ∣ξh∣≤π|\xi h| \le \pi, 存在复平面的单位圆上关于虚轴对称的两点 λ1,λ2\lambda_1, \lambda_2 满足 (1.17a), 而除此以外 (1.17b) 没有其他解 (代数基本定理), 故 λ1,λ2\lambda_1, \lambda_2 是 GG 的全部特征值. 因它们的幅值均为 1, 故有 (1.22). 充分性证毕. 证毕.

    严格来说, 对于格式 (1.13), Eq. (1.24) 及与之等价的 (1.22) 对稳定性条件 (1.21) 的充分性尚未被讨论. 我们很快会看到, ∣Cr∣=1|\text{Cr}| = 1 不能保证格式 (1.13) 的稳定性.

    格式 (1.13) 稳定振荡的现象, 可用 Fourier 分析予以解释. 在格式稳定的必要条件 (1.23) 被满足的前提下, 考虑 ξ\xi -单波 (1.14). 对于 ∣Cr sin⁡(ξh)∣<1|\text{Cr} \, \sin{(\xi h)}| < 1, 复增长矩阵 (1.16) 有两个不相等的特征值, 从而可被相似对角化. 又由 (1.15) 知, 存在可逆矩阵 P=P(ξ,τ,h)P = P(\xi, \tau, h), 使

Cjn+1=P[λ1λ2]nP−1Cj1,\begin{equation*} \tag{1.25} \boldsymbol{C}_j^{n+1} = P \begin{bmatrix*} \lambda_1 & \\ & \lambda_2 \end{bmatrix*}^n P^{-1} \boldsymbol{C}_j^{1}, \end{equation*}

其中

λ1,2=±1−Cr2sin⁡2(ξh)−i Cr sin⁡(ξh)\begin{align*} \tag{1.26} \lambda_{1,2} = \pm \sqrt{1 - \text{Cr}^2 \sin^2 (\xi h)} - \mathrm{i} \, \text{Cr} \, \sin{(\xi h)} \end{align*}

是复增长矩阵 (1.16) 的两个特征值. 另一方面, ξ\xi- 单波 Cjn\boldsymbol{C}_j^n 由 (1.1) 的方程所决定的精确行为是

C~jn+1=P[λ~λ~]nP−1C~j1,\begin{equation*} \tag{1.27} \tilde{\boldsymbol{C}}_j^{n+1} = P \begin{bmatrix*} \tilde{\lambda} & \\ & \tilde{\lambda} \end{bmatrix*}^n P^{-1} \tilde{\boldsymbol{C}}_j^{1}, \end{equation*}

其中

λ~=exp⁡(−iξuτ)=exp⁡(−i ξh Cr).\begin{equation*} \tag{1.28} \tilde{\lambda} = \exp{(-\mathrm{i} \xi u \tau)} = \exp{(-\mathrm{i} \, \xi h \, \text{Cr})}. \end{equation*}

可见, 此时格式 (1.13) 稳定, 且数值行为特征值 (1.26) 与精确行为特征值 (1.28) 的幅值相同 (均为 1), 但相位不同, 即存在相位误差. 特别地, 当 ξh Cr→0\xi h \, \text{Cr} \to 0, 有 λ~→1\tilde{\lambda} \to 1, 而 λ1,2→±1\lambda_{1,2} \to \pm 1. 此时, 数值模式 1 的特征值 λ1\lambda_1 趋于精确行为, 对应物理模式 (physical mode). 而数值模式 2 的特征值 λ2\lambda_2 远离精确行为, 对应计算模式 (computational mode), 是由数值格式引入的虚假模式. 在这种极限情形下, 计算模式虽然保持特征值幅值为 1, 从而不引起不稳定, 但每积分一个时间步都要变化一次正负号, 引起数值解稳定振荡. 一般地, 对取定的 ξh\xi h, 随着 Cr\text{Cr} 从 00 变化到 11, 数值行为特征值 (1.26) 和精确行为特征值 (1.28) 在复平面的单位圆上运动 (下图).

    对于 Cr sin⁡(ξh)=±1\text{Cr} \, \sin{(\xi h)} = \pm 1 (此时 ∣Cr∣=1|\text{Cr}| = 1 且 ∣ξh∣=π/2|\xi h| = \pi / 2), 情况将有所不同. 此时, 复增长矩阵 (1.16) 成为

G=[∓2i11],\begin{equation*} \tag{1.29} G = \begin{bmatrix} \mp 2 \mathrm{i} & 1\\ 1 & \end{bmatrix}, \end{equation*}

其特征值 (1.26)

λ1=λ2=∓i=λ~\begin{equation*} \tag{1.30} \lambda_1 = \lambda_2 = \mp \mathrm{i} = \tilde{\lambda} \end{equation*}

与精确行为的特征值 (1.28) 相同. 这特征值的几何重数 (特征子空间的维数) 2−rank(λI−G)=12 - \mathop{\text{rank}}(\lambda I - G) = 1 小于代数重数 22, 故 GG 与其 Jordan 标准型

J=[∓i1∓i]\begin{equation*} \tag{1.31} J = \begin{bmatrix} \mp \mathrm{i} & 1 \\ & \mp \mathrm{i} \end{bmatrix} \end{equation*}

相似,

G=PJP−1,\begin{equation*} \tag{1.32} G = PJP^{-1}, \end{equation*}

且 PP 可取为

P=(η1,η2),\begin{equation*} \tag{1.33} P = (\eta_1, \eta_2), \end{equation*}

其中

η1=[∓i1]′,η2=[10]′\begin{equation*} \tag{1.34} \eta_1 = \begin{bmatrix} \mp \mathrm{i} & 1 \end{bmatrix}', \quad \eta_2 = \begin{bmatrix} 1 & 0 \end{bmatrix}' \end{equation*}

是方程

(G−λI)η1=0,(G−λI)η2=η1\begin{align*} \tag{1.35a} & (G - \lambda I) \eta_1 = 0, \\ \tag{1.35b} & (G - \lambda I) \eta_2 = \eta_1 \end{align*}

的一个解. 由 (1.15), 有

Cjn+1=PJnP−1Cj1=(λ~)nCj1−n(∓i)n[∓icj0cj0]′.\begin{equation*} \tag{1.36a} \begin{aligned} \boldsymbol{C}_j^{n+1} & = P J^n P^{-1} \boldsymbol{C}_j^{1} \\ & = (\tilde{\lambda})^n \boldsymbol{C}_j^1 - n (\mp \mathrm{i})^n \begin{bmatrix} \mp \mathrm{i} c_j^0 & c^0_j \end{bmatrix}'. \end{aligned} \end{equation*}

Eq. (1.36a) 右端第一项与精确行为 (1.27) 一致, 第二项是幅度不断增长 (是 O(τ−1)\mathcal{O}(\tau^{-1}) 量级的无穷大量 as τ→0+\tau \to 0^+) 的误差项. 把 (1.36a) 写成

cjn=(1−n)(∓i)ncj0=(1−n)c~j0,\begin{equation*} \tag{1.36b} c_j^n = (1 - n) (\mp \mathrm{i})^n c_j^0 = (1 - n) \tilde{c}_j^0, \end{equation*}

看出对于使 Cr sin⁡(ξh)=±1\text{Cr} \, \sin{(\xi h)} = \pm 1 的 ξ\xi -单波, 用格式 (1.13) 的积分结果没有相位误差, 但有不断增长的幅值误差, the scheme is unstable and divergent. 这是对于周期信号而言的. 而对于 hZh \mathbb{Z} 上的平方可和的信号, 其频谱密度关于波数 ξ\xi 是在 [−π/h,π/h)[-\pi / h, \pi / h) 上连续分布, 故至多可数个频谱分量的强度不影响叠加结果, 所以认为 Cr=1\mathrm{Cr} = 1 依然保证格式 (1.13) 的稳定性. 对于有限长度的非周期信号, 格式 (1.13) 需要设置合适的初始和边界条件, 上述 Fourier 分析似乎不能严格地适用?

    综上, 格式 (1.13) 稳定的充要条件是

∣Cr∣=∣u∣τ/h<1.\begin{equation*} \tag{1.37} |\text{Cr}| = |u| \tau / h < 1. \end{equation*}

当格式稳定时, 数值解中存在虚假的计算模式, 特征值 (1.26) 与精确行为特征值 (1.28) 的差异造成数值解的稳定振荡现象.

    Eq. (1.36b) 是从频谱空间得到的结果; 下面, 我们转向物理空间, 探讨 Cr=±1\mathrm{Cr} = \pm 1 时格式 (1.13) 的稳定性. 这时 Leap-frog central diff. (LF_CD) 迭代格式 (1.3) 成为

cjn+1−cj∓1n=(−1)(cj±1n−cjn−1).\begin{equation*} \tag{1.38} c_j^{n+1} - c_{j \mp 1}^n = (-1)(c_{j \pm 1}^n - c_j^{n-1}). \end{equation*}

若首个时间步的数值积分结果精确, 例如用 Cr=±1\mathrm{Cr} = \pm 1 的迎风格式使 c0c^0 和 c1c^1 还原精确解, 则用 LF_CD 方案 (1.38) 的数值积分

cjn=cj∓n0=c0(xj∓nτu)\begin{equation*} \tag{1.39} c^n_j = c^0_{j \mp n} = c^0(x_j \mp n\tau u) \end{equation*}

就是精确解.

    对于 Cr=1\mathrm{Cr} = 1 的 LF_CD 方案 (1.38), 若在首个时间步的积分引入一个周期为四个网格距的扰动, 情况将明显不同. 定义

ejn:=cjn−c~jn\begin{equation*} \tag{1.40} e_j^n := c_j^n - \tilde{c}_j^n \end{equation*}

为绝对误差, 其中 (不妨设置 u>0u > 0)

c~jn:=c0(xj−nτu)\begin{equation*} \tag{1.41} \tilde{c}_j^n := c^0(x_j - n \tau u) \end{equation*}

是一维平流问题 (1.1) 的精确解. 由 (1.38) 得误差传递关系

ejn+1−ej−1n=(−1)(ej+1n−ejn−1)=(−1)n(ej+n1−ej+n−10),\begin{align*} \tag{1.42a} e_j^{n+1} - e_{j-1}^n & = (-1)(e^n_{j+1} - e^{n-1}_j) \\ \tag{1.42b} & = (-1)^n (e^1_{j+n} - e^0_{j+n-1}), \end{align*}

进而

ejn+1=ej−n−10+∑k=0n(−1)kej−n+2k1.\begin{equation*} \tag{1.42c} e_j^{n+1} = e^0_{j-n-1} + \sum_{k = 0}^{n} {(-1)^k e^1_{j-n+2k}}. \end{equation*}

设初始状态被精确离散 (i.e., e0=0e^0 = 0), 且首个时间步的积分引入误差

ej1=Acos⁡(jπ/2)+Bsin⁡(jπ/2),\begin{equation*} \tag{1.43} e_j^1 = A \cos{(j \pi / 2)} + B \sin{(j \pi / 2)}, \end{equation*}

则由 (1.42c)(1.43) 得

ejn+1=∑k=0n(−1)k[Acos⁡(kπ+(j−n)π/2)+Bsin⁡(kπ+(j−n)π/2)]=∑k=0n[Acos⁡((j−n)π/2)+Bsin⁡((j−n)π/2)],\begin{equation*} \tag{1.44} \begin{aligned} e_j^{n+1} & = \sum_{k = 0}^{n} { (-1)^k \left[ A \cos{(k\pi + (j-n) \pi / 2)} + B \sin{(k\pi + (j-n) \pi / 2)} \right] } \\ & = \sum_{k = 0}^{n} { \left[ A \cos{((j-n) \pi / 2)} + B \sin{((j-n) \pi / 2)} \right] }, \end{aligned} \end{equation*}

进而长度为 MM 的误差向量 ene^n 的 L2L_2-norm 满足

M−1/2∥en∥2→n(A2+B2)/2,(M→+∞),\begin{equation*} \tag{1.45} M^{-1/2} \| e^n \|_2 \to n \sqrt{(A^2 + B^2) / 2}, \quad (M \to +\infty), \end{equation*}

即误差按 O(n)=O(τ−1)\mathcal{O}(n) = \mathcal{O}(\tau^{-1}) 量级增长, LF_CD scheme (1.13) is un-stable and in-consistent. 而若首个时间步积分没有误差, 则用 Cr=1\mathrm{Cr} = 1 的 LF_CD scheme (1.13) 将得到精确解 (1.39); 但这种临界稳定状态很危险, 首步积分稍有误差就可能 blows up, 所以最好让更严格的 (1.22) 即 (1.37) 成立以保证稳定性.

    从时域分析得到的 Eq. (1.45), 与从 Fourier 频域分析得到 Eq. (1.39) 意义相近, 都是说 Cr sin⁡(ξh)=±1\mathrm{Cr} \, \sin{(\xi h)} = \pm 1 的单色波经过 LF_CD 格式 (1.3) 会 blows up. 然而, 时域和频域结果在细节上似乎有些差异. 例如, 时域结果 (1.39)(1.45) 能解释本文第 1(a) 节图 (c) 为何在临界稳定的情况下能得到精确解, 而这似乎是频域分析结果 Eq. (1.39) 所不能直接做到的. 对这种差异的数学和物理解释, 可能是重要的, which 关系对 von Neumann 分析方法的准确理解. 我在此呼吁进一步讨论, 请求有兴趣的读者提供帮助! 文末附有联系方式. 下面的简单试验也许可作为讨论的开始.

    为验证 (1.36) 或 (1.45), 设计了下述实验. 对 LF_CD 格式 (1.13), 提供零初值, 采取周期边界, 周期是扰动信号周期的 22+k, k∈N∗2^{2 + k}, \, k \in \mathbb{N}^* 倍. 首步积分采取迎风格式 (参见 Hw 2 - Problem 2), 对首步积分结果添加扰动信号 (1.43), 取不同的 Cr\text{Cr} 作积分. 将扰动改为 (1.43) 的倍频或半频信号, 重复实验. 结果示于下图. 试验结果表明, LF_CD 格式 (1.13) 当取 Cr=1\text{Cr} = 1 时, 对于 ξh=π/2\xi h = \pi / 2 的中频扰动信号不稳定, 误差幅值与积分时间成正比; 对于 ξh=π/4\xi h = \pi / 4 的低频信号和 ξh=π/1\xi h = \pi / 1 的高频信号稳定 (下图 (c), 下下图 (c) 中误差很小, 但也表现出正比于积分时间的增长趋势, 这如何解释?). 这就验证了 (1.45). 但如何理解 (1.36)?

fig_disturb_test_4h

fig_disturb_test_2h

fig_disturb_test_8h

    (d) 为缓解 Leap-frog central diff. 格式 (1.13) 的稳定振荡 (computational mode) 现象, 一种方案是引入 Robert–Asselin (RA) time filter (Robert[2], 1966; Asselin[3], 1972). 具体来说, 就是在每次完成 tn+1t_{n+1} 时步的积分之后, 执行赋值语句 (Ferziger[4] et al., 2020, Eq. (6.47))

cin←cin+ν2(cin+1−2cin+cin−1)=cin+df.\begin{equation*} \tag{1.46} c^n_i \leftarrow c_i^n + \frac{\nu}{2}\left( c_{i}^{n+1} - 2c_i^{n} + c_{i}^{n-1} \right) = c_i^n + d_f. \end{equation*}

针对 RA time filter (1.46) 引入的 undesired numerical damping 和精度下降的问题, Williams[5] (2009) 对 (1.38) 提出改进 (Ferziger[4:1] et al., 2020, Eq. (6.48)(6.49)),

cin←cin+αdf,cin+1←cin+1−(1−α)df.\begin{align*} \tag{1.47a} & c^n_i \leftarrow c_i^n + \alpha d_f, \\ \tag{1.47b} & c^{n+1}_i \leftarrow c_i^{n+1} - (1 - \alpha) d_f. \end{align*}

通常, 滤波器参数 ν\nu 取 0.01–0.2, 且取 α≥0.5\alpha \ge 0.5 以保证稳定性. Eq. (1.46)(1.47) 的几何意义示于下图 (Williams, 2009, his Fig. 1).

Fig. 1. of Williams (2009, doi: 10.1175/2009MWR2724.1)

    用 Matlab 编程, 为格式 (1.3) 添加 RA time filter (1.46). 数值积分结果表明, 当 ν=0.05\nu = 0.05, 虚假振荡现象略有改善, 但临界稳定 Cr\mathrm{Cr} 由 11 左移至约 0.980.98. 若增大 ν\nu 至 ν=1\nu = 1, 则虚假振荡被进一步消除, 但物理信号也被过度衰减, 且临界稳定 Cr\mathrm{Cr} 进一步左移至约 0.60, 使时间步长的取值被限制. 若引入 Williams (2009) 的改进 (1.47), 则物理信号被过度衰减的问题可获改善, 但虚假振荡的消除效果亦被削弱. 可通过优化滤波器参数 ν\nu 和 α\alpha 的选取, 使数值积分取得更优效果. 该实验体现了多参数联合调优在数值积分性能优化工作中的意义.

RA time filter with ν = 0.05 and α = 1

RA time filter with ν = 1 and α = 1

RA time filter with ν = 0.53 and α = 1

RA time filter with ν = 0.2 and α = 53

Problem 2

    Consider the following numerical method which is proposed for solving ∂u∂t=c∂u∂x\frac{\partial u}{\partial t} = c \frac{\partial u}{\partial x},

1Δt[ujn+1−12(uj+1n+uj−1n)]=c2Δx(uj+1n−uj−1n).\frac{1}{\Delta t} \left[ u_j^{n+1} - \frac{1}{2} \left(u^n_{j+1} + u^n_{j-1} \right) \right] = \frac{c}{2 \Delta x} \left( u^n_{j+1} - u^n_{j-1} \right).

Use von Neumann stability analysis to find the range of Courant numbers for which this method is stable. Perform modified equation analysis on this method. Is the method consistent, i.e. does it reduce to the original PDE as Δx\Delta x and Δt\Delta t each go to zero?

    满足格式

1τ[ujn+1−12(uj+1n+uj−1n)]=c2h(uj+1n−uj−1n)\begin{equation*} \tag{2.1} \frac{1}{\tau} \left[ u_j^{n+1} - \frac{1}{2} \left(u^n_{j+1} + u^n_{j-1} \right) \right] = \frac{c}{2 h} \left( u^n_{j+1} - u^n_{j-1} \right) \end{equation*}

的 u~∈C6(R)\tilde{u} \in C^6(\mathbb{R}) 成立 modified equation

u~t−cu~x=τ2u~t2−h22τu~x2+h26cu~x3+O(τ2+h4τ−1+h4),\begin{equation*} \tag{2.2} \tilde{u}_t - c\tilde{u}_x = \frac{\tau}{2} \tilde{u}_{t^2} - \frac{h^2}{2 \tau} \tilde{u}_{x^2} + \frac{h^2}{6} c \tilde{u}_{x^3} + \mathcal{O}(\tau^2 + h^4 \tau^{-1} + h^4), \end{equation*}

故格式 (2.1) is consistent iff. (τ,h)(\tau, h) 沿 h2=o(τ)  as  τ→0h^2 = o(\tau) \; \mathrm{as} \; \tau \to 0 的路径趋于 (0,0)(0,0).

    格式 (2.1) 的复增长因子

G=cos⁡(ξh)+icτhsin⁡(ξh)\begin{equation*} \tag{2.3} G = \cos{(\xi h)} + \mathrm{i} \frac{c \tau}{h} \sin{(\xi h)} \end{equation*}

满足加强的 von Neumann 条件 (1.22)

∀∣ξh∣≤π:  ∣G∣2=1+((cτ/h)2−1)sin⁡2(ξh)≤1\begin{equation*} \tag{2.4} \forall |\xi h| \le \pi : \; |G|^2 = 1 + ((c \tau / h)^2 - 1) \sin^2{(\xi h)} \le 1 \end{equation*}

iff. Courant number 的平方

Cr2:=(cτ/h)2≤1.\begin{equation*} \tag{2.5} \mathrm{Cr}^2 := (c \tau / h)^2 \le 1. \end{equation*}

Acknowledgement

    I am grateful to my fellow PhD student, Dan Cao (曹丹), whose insightful discussions significantly enriched this research.

Contact Information

Appendix

Matlab code for problem 1

hw3_1.m 

%% hw3_1.m
% Description: Solving the 1-D advection problem: c_t + u c_x = 0
% Author: Guorui Wei (危国锐) (313017602@qq.com)
% Created at: Nov. 7, 2024
% Last modified: Nov. 9, 2024
%

clear; clc; close all

addpath("..\hw2\")
if ~isfolder(".\fig\")
    mkdir(".\fig\")
end

%% Problem 1(a)

Cr_list = [.50, .99, 1.00, 1.01]; % Cr = u delta_t / delta_x
[f_list, x_grid, t_list] = hw3_1_unit("time_leap_frog_RA_space_central_diff", Cr_list, true, nu=0, alpha=1);

%% Problem 1(d)

Cr_list = [.50, .75, .975, 1]; % Cr = u delta_t / delta_x
[f_list, x_grid, t_list] = hw3_1_unit("time_leap_frog_RA_space_central_diff", Cr_list, true, nu=.05, alpha=1);

Cr_list = [.50, .75, .975, 1]; % Cr = u delta_t / delta_x
[f_list, x_grid, t_list] = hw3_1_unit("time_leap_frog_RA_space_central_diff", Cr_list, true, nu=.05, alpha=.53);

Cr_list = [.50, .25, .57, .60]; % Cr = u delta_t / delta_x
[f_list, x_grid, t_list] = hw3_1_unit("time_leap_frog_RA_space_central_diff", Cr_list, true, nu=1, alpha=1);

Cr_list = [.50, .25, .57, .60]; % Cr = u delta_t / delta_x
[f_list, x_grid, t_list] = hw3_1_unit("time_leap_frog_RA_space_central_diff", Cr_list, true, nu=1, alpha=.53);

Cr_list = [.50, .25, .57, .60]; % Cr = u delta_t / delta_x
[f_list, x_grid, t_list] = hw3_1_unit("time_leap_frog_RA_space_central_diff", Cr_list, true, nu=.2, alpha=.53);

%% 首步扰动试验

clc;clear;close all;

scheme_name = "time_leap_frog_RA_space_central_diff";
Cr_list = [.5, .99, 1, 1.01];
delta_x = 0.25;
x_range = {[0, delta_x*4*(2^2)], };
% init_func = @(x) double(x > 0 & ~(x > .5));
init_func = @(x) zeros(size(x));
bndry_func = @(x, t) "periodic";
t_query = [0, 1, 2, 4, 8, 16];
ra_opts.nu = 0;
ra_opts.alpha = 1;
flag_save = true;

%%% 周期 4h

disturb_func_first_step = @(x) 1*(sin(2*pi / 4 / delta_x * x) + cos(2*pi / 4 / delta_x * x));
disturb_func_disp_str = "$f_{\rm{dist}} = 1 \cdot (\sin{(2 \pi x / (4h))} + \cos{(2 \pi x / (4h))})$";
fig_name = "fig_disturb_test_4h";

[f_list, x_grid, t_list] = disturb_test_unit(scheme_name, Cr_list, x_range, init_func, bndry_func, disturb_func_first_step, disturb_func_disp_str, delta_x, t_query, flag_save, fig_name, nu=ra_opts.nu, alpha=ra_opts.alpha);

%%% 周期 2h (Nyquist)

disturb_func_first_step = @(x) 1*(sin(2*pi / 2 / delta_x * x) + cos(2*pi / 2 / delta_x * x));
disturb_func_disp_str = "$f_{\rm{dist}} = 1 \cdot (\sin{(2 \pi x / (2h))} + \cos{(2 \pi x / (2h))})$";
fig_name = "fig_disturb_test_2h";

[f_list, x_grid, t_list] = disturb_test_unit(scheme_name, Cr_list, x_range, init_func, bndry_func, disturb_func_first_step, disturb_func_disp_str, delta_x, t_query, flag_save, fig_name, nu=ra_opts.nu, alpha=ra_opts.alpha);

%%% 周期 8h

disturb_func_first_step = @(x) 1*(sin(2*pi / 8 / delta_x * x) + cos(2*pi / 8 / delta_x * x));
disturb_func_disp_str = "$f_{\rm{dist}} = 1 \cdot (\sin{(2 \pi x / (8h))} + \cos{(2 \pi x / (8h))})$";
fig_name = "fig_disturb_test_8h";

[f_list, x_grid, t_list] = disturb_test_unit(scheme_name, Cr_list, x_range, init_func, bndry_func, disturb_func_first_step, disturb_func_disp_str, delta_x, t_query, flag_save, fig_name, nu=ra_opts.nu, alpha=ra_opts.alpha);

%% Function: Unit

function [f_list, x_grid, t_list] = hw3_1_unit(scheme_name, Cr_list, flag_save, ra_opts)
    arguments (Input)
        scheme_name = "time_leap_frog_RA_space_central_diff";
        Cr_list = [.5, .99, 1, 1.01]; 
        flag_save = false;
        ra_opts.nu = 0;
        ra_opts.alpha = 1;
    end
    
    solver = AdvecProb2D();

    %%% Solver params (1-D Advection problem)

    velocity = [1, ];
    t_start = 0;
    x_range = {[0, 4], };
    init_func = @(x) double(x > 0 & ~(x > .5));
    bndry_func = @(x, t) double(x > x_range{1}(1) & x < x_range{1}(2));
    delta_x = 0.05;
    t_query = [0, 1, 2, 3];

    %%% Solving

    % Cr_list = [.5, .99, 1, 1.01]; % Cr = u delta_t / delta_x
    f_list = cell(length(Cr_list), 1);
    x_grid = f_list;
    t_list = x_grid;
    for test_idx = 1:length(Cr_list)
        Cr = Cr_list(test_idx);
        delta_t =  Cr * delta_x / velocity(1);
        params_struct = AdvecProb2D.prepare_params(velocity, t_start, x_range, init_func, bndry_func, delta_t, delta_x, t_query);
        solver.reset(params_struct);
        solver.solve(scheme_name, nu=ra_opts.nu, alpha=ra_opts.alpha);
        [f_list{test_idx}, x_grid{test_idx}, t_list{test_idx}, ~] = solver.get_solution();
    end

    %%% Create figure
    switch scheme_name
        case "time_forward_Euler_space_upwind"
            fig_name_prefix = "fig_upwind";
            scheme_name_short = "upwind";
        case "time_forward_Euler_space_central_diff"
            fig_name_prefix = "fig_centdif";
            scheme_name_short = "central diff.";
        case "time_leap_frog_RA_space_central_diff"
            fig_name_prefix = sprintf("fig_LF_RA_CD_nu_%.2g_alpha_%.2g", ra_opts.nu, ra_opts.alpha);
            scheme_name_short = "LF\_RA\_CD";
        otherwise
            fig_name_prefix = "fig_undefined";
            scheme_name_short = "unknown scheme";
    end

    t_fig = figure(Name=sprintf("%s_Cr_compare_1", fig_name_prefix));

    % set figure size
    UNIT_ORIGINAL = t_fig.Units;
    t_fig.Units = "centimeters";
    t_fig.Position = [3, 3, 16, 16];
    t_fig.Units = UNIT_ORIGINAL;

    % plot
    t_TCL = tiledlayout(t_fig, floor(length(Cr_list) / 2), 2, TileSpacing="compact", Padding="compact");
    xlabel(t_TCL, "$x$", Interpreter="latex", FontSize=10.5);
    ylabel(t_TCL, "(numerical, " + scheme_name_short + ") $f$ (subjected to $f_{t} + u f_x = 0, \; f|_{t = 0} = f^0, \; f|_{x \in \partial \Omega} = 0$)", Interpreter="latex", FontSize=10.5);
    if scheme_name == "time_leap_frog_RA_space_central_diff"
        title(t_TCL, {"$c^n_i \leftarrow c_i^n + \alpha d_f, \; c^{n+1}_i \leftarrow c_i^{n+1} - (1 - \alpha) d_f, \; d_f = \frac{\nu}{2} \left( c_{i}^{n+1} - 2c_i^{n} + c_{i}^{n-1} \right)$"; sprintf("$\\nu$ = %.3g, $\\alpha$ = %.3g", ra_opts.nu, ra_opts.alpha)}, Interpreter="latex", FontSize=10.5);
    end

    for test_idx = 1:length(Cr_list)
        t_axes = nexttile(t_TCL, test_idx);
        hold(t_axes, "on");
        for t_query_idx = 1:length(t_list{test_idx})
            plot(t_axes, x_grid{test_idx}{1}, f_list{test_idx}{t_query_idx}, LineWidth=1.0, DisplayName=sprintf("$t = %.2f$", t_list{test_idx}(t_query_idx)))
        end
        set(t_axes, FontName="Times New Roman", FontSize=10.5, Box="on", TickLabelInterpreter="latex", ...
            Tag=sprintf("(%c) $\\mathrm{Cr} = %g$", char('a' - 1 + test_idx), Cr_list(test_idx)));
        legend(t_axes, Location="northwest", Interpreter="latex", FontSize=10.5, Box="off");
    end

    for t_axes = findobj(t_TCL, 'Type', "Axes", {'-regexp', 'Tag', "^\([a-z]+\)"}).'
        str_{1} = "\bf " + t_axes.Tag;
        t_txt_box = annotation(t_fig, "textbox", String=str_, Position=[t_axes.Position([1, 2]) + t_axes.Position([3, 4]), .1, .1], FontSize=10.5, Interpreter="latex", LineStyle="none", HorizontalAlignment="right", VerticalAlignment="top");
        UNIT_ORIGINAL = t_txt_box.Units;
        t_txt_box.Units = "points";
        t_txt_box.Position = [t_txt_box.Position([1,2]) - [10.5*10, 10.5*1.5], 10.5*10, 10.5*1.5];
        t_txt_box.Units = UNIT_ORIGINAL;
    end

    if flag_save
        print(t_fig, ".\fig\" + t_fig.Name + ".svg", "-vector", "-dsvg")
    end

    %%% Fig. Cr compare fixed t
    t_fixed = 2;
    t_fig = figure(Name=sprintf("%s_Cr_compare_2", fig_name_prefix));

    % set figure size
    UNIT_ORIGINAL = t_fig.Units;
    t_fig.Units = "centimeters";
    t_fig.Position = [3, 3, 12, 12];
    t_fig.Units = UNIT_ORIGINAL;

    % plot
    t_TCL = tiledlayout(t_fig, 1, 1, TileSpacing="compact", Padding="compact");
    if scheme_name == "time_leap_frog_RA_space_central_diff"
        title(t_TCL, {"$c^n_i \leftarrow c_i^n + \alpha d_f, \; c^{n+1}_i \leftarrow c_i^{n+1} - (1 - \alpha) d_f, \; d_f = \frac{\nu}{2} \left( c_{i}^{n+1} - 2c_i^{n} + c_{i}^{n-1} \right)$"; sprintf("$\\nu$ = %.3g, $\\alpha$ = %.3g", ra_opts.nu, ra_opts.alpha)}, Interpreter="latex", FontSize=10.5);
    end

    t_axes = nexttile(t_TCL, 1);

    test_idx = find(Cr_list == 1, 1);
    if ~isempty(test_idx)
        [~, t_query_idx] = min(t_list{test_idx} - t_fixed, [], ComparisonMethod="abs");
        plot(t_axes, x_grid{test_idx}{1}, f_list{test_idx}{t_query_idx}, "-k", LineWidth=1.0, DisplayName=sprintf("$\\mathrm{Cr} = %g, \\; t = %.2f$", Cr_list(test_idx), t_list{test_idx}(t_query_idx)));
    end
    hold(t_axes, "on");
    for test_idx = 1:length(Cr_list)
        if Cr_list(test_idx) == 1
            continue
        end
        [~, t_query_idx] = min(t_list{test_idx} - t_fixed, [], ComparisonMethod="abs");
        plot(t_axes, x_grid{test_idx}{1}, f_list{test_idx}{t_query_idx}, LineWidth=1.0, DisplayName=sprintf("$\\mathrm{Cr} = %g, \\; t = %.2f$", Cr_list(test_idx), t_list{test_idx}(t_query_idx)));
    end

    set(t_axes, FontName="Times New Roman", FontSize=10.5, Box="on", TickLabelInterpreter="latex", XLimitMethod="tight");
    legend(t_axes, Location="best", Interpreter="latex", FontSize=10.5, Box="off");

    xlabel(t_axes, "$x$", Interpreter="latex", FontSize=10.5);
    ylabel(t_axes, "(numerical, " + scheme_name_short + ") $f$", Interpreter="latex", FontSize=10.5);

    if flag_save
        print(t_fig, ".\fig\" + t_fig.Name + ".svg", "-vector", "-dsvg")
    end

end

%% Func: 首步扰动试验

function [f_list, x_grid, t_list] = disturb_test_unit(scheme_name, Cr_list, x_range, init_func, bndry_func, disturb_func_first_step, disturb_func_disp_str, delta_x, t_query, flag_save, fig_name, ra_opts)
    arguments (Input)
        scheme_name = "time_leap_frog_RA_space_central_diff";
        Cr_list = [.5, .99, 1, 1.01];
        x_range = {[0, 2], };
        init_func = @(x) 5*double(x > 0 & ~(x > .5));
        bndry_func = @(x, t) "periodic";
        disturb_func_first_step = @(x) sin(pi / 2 / obj.delta_x * x) + cos(pi / 2 / obj.delta_x * x);
        disturb_func_disp_str = "$f_{\rm{dist}} = \sin{(\pi x / (2h))} + \cos{(\pi x / (2h))}$";
        delta_x = 0.05;
        t_query = [0, 1, 2, 3];
        flag_save = false;
        fig_name = "fig_disturb_test";
        ra_opts.nu = 0;
        ra_opts.alpha = 1;
    end
    
    solver = AdvecProb2D();

    %%% Solver params (1-D Advection problem)

    velocity = [1, ];
    t_start = 0;
    % x_range = {[0, 4], };
    % init_func = @(x) double(x > 0 & ~(x > .5));
    % bndry_func = @(x, t) double(x > x_range{1}(1) & x < x_range{1}(2));
    % delta_x = 0.05;
    % t_query = [0, 1, 2, 3];

    %%% Solving

    % Cr_list = [.5, .99, 1, 1.01]; % Cr = u delta_t / delta_x
    f_list = cell(length(Cr_list), 1);
    x_grid = f_list;
    t_list = x_grid;
    for test_idx = 1:length(Cr_list)
        Cr = Cr_list(test_idx);
        delta_t =  Cr * delta_x / velocity(1);
        params_struct = AdvecProb2D.prepare_params(velocity, t_start, x_range, init_func, bndry_func, delta_t, delta_x, t_query);
        solver.reset(params_struct);
        solver.solve(scheme_name, disturb_func_first_step, nu=ra_opts.nu, alpha=ra_opts.alpha);
        [f_list{test_idx}, x_grid{test_idx}, t_list{test_idx}, ~] = solver.get_solution();
    end

    %%% Create figure
    switch scheme_name
        case "time_forward_Euler_space_upwind"
            scheme_name_short = "upwind";
        case "time_forward_Euler_space_central_diff"
            scheme_name_short = "central diff.";
        case "time_leap_frog_RA_space_central_diff"
            scheme_name_short = "LF\_RA\_CD";
        otherwise
            scheme_name_short = "unknown scheme";
    end

    t_fig = figure(Name=fig_name);

    % set figure size
    UNIT_ORIGINAL = t_fig.Units;
    t_fig.Units = "centimeters";
    t_fig.Position = [3, 3, 16, 16];
    t_fig.Units = UNIT_ORIGINAL;

    % plot
    t_TCL = tiledlayout(t_fig, floor(length(Cr_list) / 2), 2, TileSpacing="compact", Padding="compact");
    xlabel(t_TCL, "$x$", Interpreter="latex", FontSize=10.5);
    ylabel(t_TCL, "(numerical, " + scheme_name_short + ") $f$ (subjected to $f_{t} + u f_x = 0, \; f|_{t = 0} = f^0, \; f|_{x \in \partial \Omega} = 0$)", Interpreter="latex", FontSize=10.5);
    if scheme_name == "time_leap_frog_RA_space_central_diff"
        title(t_TCL, {"$c^n_i \leftarrow c_i^n + \alpha d_f, \; c^{n+1}_i \leftarrow c_i^{n+1} - (1 - \alpha) d_f, \; d_f = \frac{\nu}{2} \left( c_{i}^{n+1} - 2c_i^{n} + c_{i}^{n-1} \right)$"; ...
                      sprintf("$\\nu$ = %.3g, $\\alpha$ = %.3g, $\\quad$ %s", ra_opts.nu, ra_opts.alpha, disturb_func_disp_str) ...
                      }, Interpreter="latex", FontSize=10.5);
    end

    for test_idx = 1:length(Cr_list)
        t_axes = nexttile(t_TCL, test_idx);
        hold(t_axes, "on");
        for t_query_idx = 1:length(t_list{test_idx})
            plot(t_axes, x_grid{test_idx}{1}, f_list{test_idx}{t_query_idx}, LineWidth=1.0, DisplayName=sprintf("$t = %.2f$", t_list{test_idx}(t_query_idx)))
        end
        set(t_axes, FontName="Times New Roman", FontSize=10.5, Box="on", TickLabelInterpreter="latex", ...
            Tag=sprintf("(%c) $\\mathrm{Cr} = %g$", char('a' - 1 + test_idx), Cr_list(test_idx)));
        legend(t_axes, Location="northwest", Interpreter="latex", FontSize=10.5, Box="off");
    end

    for t_axes = findobj(t_TCL, 'Type', "Axes", {'-regexp', 'Tag', "^\([a-z]+\)"}).'
        str_{1} = "\bf " + t_axes.Tag;
        t_txt_box = annotation(t_fig, "textbox", String=str_, Position=[t_axes.Position([1, 2]) + t_axes.Position([3, 4]), .1, .1], FontSize=10.5, Interpreter="latex", LineStyle="none", HorizontalAlignment="right", VerticalAlignment="top");
        UNIT_ORIGINAL = t_txt_box.Units;
        t_txt_box.Units = "points";
        t_txt_box.Position = [t_txt_box.Position([1,2]) - [10.5*10, 10.5*1.5], 10.5*10, 10.5*1.5];
        t_txt_box.Units = UNIT_ORIGINAL;
    end

    if flag_save
        print(t_fig, ".\fig\" + t_fig.Name + ".svg", "-vector", "-dsvg")
    end
end

AdvecProb2D.m 

%% AdvecProb2D.m
% Description: Solving the 1-D or 2-D advection problem: c_t + u c_x + v c_y = 0, u, v constants
% Author: Guorui Wei (危国锐) (313017602@qq.com)
% Created at: Oct. 21, 2024
% Last modified: Nov. 9, 2024
%

%%

classdef AdvecProb2D < handle
    properties (Constant)
        PARAMS_REQUIRED = ["velocity","t_start","x_range","init_func","bndry_func", "delta_t", "delta_x", "t_query"];
        SCHEME_NAME = ["time_forward_Euler_space_upwind", "time_forward_Euler_space_central_diff", "time_forward_Euler_space_CTU", "time_leap_frog_RA_space_central_diff"];
    end

    properties (Access=private)
        % solver params
        velocity
        t_start
        x_range
        init_func
        bndry_func
        delta_t
        delta_x
        t_query

        % results
        x_grid
        t_list
        f_list
        scheme_name
    end

    methods (Static, Access=public)
        function params_struct = prepare_params(velocity, t_start, x_range, init_func, bndry_func, delta_t, delta_x, t_query)
            arguments (Input)
                velocity = [1, 1];
                t_start = 0;
                x_range = {[0, 1], [0, 1]};
                init_func = @(x, y) (((x - 1/2).^2 + (y - 1/2).^2)*16 < 1) .* (1 + cos(pi * 4*sqrt( (x - 1/2).^2 + (y - 1/2).^2 ))) / 2;
                bndry_func = @(x, y, t) "periodic";
                delta_t = .0125 / 5;
                delta_x = [.0250, .0250] / 5;
                t_query = [0, .5, 1] * 5;  
            end

            params_struct.velocity = velocity;
            params_struct.t_start = t_start;
            params_struct.x_range = cellfun(@(list_) sort(list_, "ascend"), x_range, UniformOutput=false); 
            params_struct.init_func = init_func;
            params_struct.bndry_func = bndry_func;
            params_struct.delta_t = delta_t;
            params_struct.delta_x = delta_x;
            params_struct.t_query = sort(t_query(t_query >= t_start), "ascend");
        end
    end

    methods (Access=public)
        function obj = AdvecProb2D(varargin)
            if nargin < 1
                obj.reset(AdvecProb2D.prepare_params());
                return
            end
            obj.reset(varargin{:});
        end

        function obj = reset(obj, params_struct)
            field_missing_flag = ~isfield(params_struct, AdvecProb2D.PARAMS_REQUIRED);
            if any(field_missing_flag)
                warning("Invalid params_struct. Params not changed. Missing field: %s", join(AdvecProb2D.PARAMS_REQUIRED(field_missing_flag), ", "));
                return
            end

            params_struct = AdvecProb2D.prepare_params(params_struct.velocity, params_struct.t_start, params_struct.x_range, params_struct.init_func, params_struct.bndry_func, params_struct.delta_t, params_struct.delta_x, params_struct.t_query);
            for field_name = AdvecProb2D.PARAMS_REQUIRED
                obj.(field_name) = params_struct.(field_name);
            end

            obj.t_list = nan(size(obj.t_query));
            obj.f_list = cell(size(obj.t_query));
        end
    
        function [f_list, x_grid, t_list, scheme_name] = get_solution(obj)
            f_list = obj.f_list;
            x_grid = obj.x_grid;
            t_list = obj.t_list;
            scheme_name = obj.scheme_name;
        end

        function obj = solve(obj, scheme_name, disturb_func_first_step, ra_opts)
            arguments (Input)
                obj 
                scheme_name
                disturb_func_first_step = @(x) 0;
                ra_opts.nu = 0;
                ra_opts.alpha = 1;
            end

            if isscalar(obj.velocity)
                [obj.f_list, obj.x_grid, obj.t_list] = obj.solve_1D(scheme_name, disturb_func_first_step, nu=ra_opts.nu, alpha=ra_opts.alpha);
                return
            end
            if length(obj.velocity) == 2
                [obj.f_list, obj.x_grid, obj.t_list] = obj.solve_2D(scheme_name);
                return
            end
        end
    end

    methods (Access=private)
        function [f_list, x_grid, t_list] = solve_1D(obj, scheme_name, disturb_func_first_step, ra_opts)
            arguments (Input)
                obj
                scheme_name;
                disturb_func_first_step = @(x) sin(pi / 2 / obj.delta_x * x) + cos(pi / 2 / obj.delta_x * x);
                ra_opts.nu = 0; % Robert-Asselin filter params
                ra_opts.alpha = 1;
            end

            f_list = cell(size(obj.t_query));
            t_list = nan(size(obj.t_query));
            x_grid = cellfun(@(range_1D, delta_1D) (range_1D(1): delta_1D: range_1D(2)).', obj.x_range, num2cell(obj.delta_x), UniformOutput=false);

            %%% scheme I (2-level)

            if ismember(scheme_name, ["time_forward_Euler_space_upwind", "time_forward_Euler_space_central_diff"])
                obj.scheme_name = scheme_name;

                Cr = abs(obj.velocity) * obj.delta_t ./ obj.delta_x;
                f_last = obj.init_func(x_grid{1});
                t_last = obj.t_start;
                t_now = t_last;
                t_query_next_idx = find(obj.t_query > obj.t_start, 1, "first");

                % 若请求时刻不晚于初始时刻, 则无需计算, 直接返回初值
                if isempty(t_query_next_idx)
                    t_list(1) = t_last;
                    f_list{1} = f_last;

                    fprintf("%s\n\tscheme_name = %s,\n\tFinished.\n\tt_list = [%s]\n", datetime('now'), scheme_name, join(string(t_list), ", "));
                    return
                end
                
                % 用于记录进度
                REC_TIME_INT = 5; % seconds
                time_start = tic;
                t_now_last_rec = t_last;

                flag_first_step = true;
                while t_now < obj.t_query(end)
                    t_now = t_last + obj.delta_t;

                    % 输出进度
                    
                    if toc(time_start) > REC_TIME_INT
                        time_remained = REC_TIME_INT * (obj.t_query(end) - t_now) / (t_now - t_now_last_rec);
                        fprintf("%s\n\tscheme_name = %s,\n\ttime_remained = %.2e,\n\tspace_dim = (%s), Cr = (%s),\n\tt_now = %.2e, t_query_next = %.2e, t_query_end = %.2e\n", datetime('now'), scheme_name, time_remained, join(string(size(X)), ", "), join(string(Cr), ", "), t_now, obj.t_query(t_query_next_idx), obj.t_query(end));
                        time_start = tic;
                        t_now_last_rec = t_now;
                    end

                    % add halo points, advance one time step
                    if scheme_name == "time_forward_Euler_space_upwind" && obj.velocity(1) > 0
                        if string(obj.bndry_func(0, 0)) == "periodic"
                            halo_left = f_last(end-1);
                        else
                            halo_left = obj.bndry_func(x_grid{1}(1) - obj.delta_x(1), t_last);
                        end
                        halo_right = [];
                        f_last_expanded = [halo_left; f_last; halo_right];
                        f_now = (1 - Cr(1))*f_last_expanded(1 + length(halo_left) : end-length(halo_right)) + Cr(1)*f_last_expanded(length(halo_left) : -1 + end - length(halo_right));
                    elseif scheme_name == "time_forward_Euler_space_upwind"
                        halo_left = [];
                        if string(obj.bndry_func(0, 0)) == "periodic"
                            halo_right = f_last(2);
                        else
                            halo_right = obj.bndry_func(x_grid{1}(end) + obj.delta_x(1), t_last);
                        end
                        f_last_expanded = [halo_left; f_last; halo_right];
                        f_now = (1 - Cr(1))*f_last_expanded(1 + length(halo_left) : end-length(halo_right)) + Cr(1)*f_last_expanded(2 + length(halo_left) : 1 + end-length(halo_right));
                    elseif scheme_name == "time_forward_Euler_space_central_diff"
                        if string(obj.bndry_func(0, 0)) == "periodic"
                            halo_left = f_last(end-1);
                            halo_right = f_last(2);
                        else
                            halo_left = obj.bndry_func(x_grid{1}(1) - obj.delta_x(1), t_last);
                            halo_right = obj.bndry_func(x_grid{1}(end) + obj.delta_x(1), t_last);
                        end
                        f_last_expanded = [halo_left; f_last; halo_right];
                        f_now = f_last_expanded(1 + length(halo_left) : end-length(halo_right)) + sign(obj.velocity)*Cr(1)/2*(f_last_expanded(length(halo_left) : -1 + end-length(halo_right)) - f_last_expanded(2 + length(halo_left) : 1 + end-length(halo_right)));
                    end

                    % 添加首步扰动
                    if flag_first_step
                        flag_first_step = false;
                        f_now = f_now + disturb_func_first_step(x_grid{1});
                    end

                    % store solution at queried time.
                    if ~((obj.t_query(t_query_next_idx) - t_last)*(obj.t_query(t_query_next_idx) - t_now) > 0)
                        if abs(obj.t_query(t_query_next_idx) - t_now) < abs(obj.t_query(t_query_next_idx) - t_last)
                            t_list(t_query_next_idx) = t_now;
                            f_list{t_query_next_idx} = f_now;
                        else
                            t_list(t_query_next_idx) = t_last;
                            f_list{t_query_next_idx} = f_last;
                        end
                        t_query_next_idx = t_query_next_idx + 1;
                    end

                    % prepare for next time step
                    t_last = t_now;
                    f_last = f_now;
                end % end of `while t_now < obj.t_query(end)`

                fprintf("%s\n\tscheme_name = %s,\n\tFinished.\n\tt_list = [%s]\n", datetime('now'), scheme_name, join(string(t_list), ", "));
                return
            end % end of `if ismember(scheme_name, ["time_forward_Euler_space_upwind", "time_forward_Euler_space_central_diff"])`

            %%% scheme II (3-level)

            if ismember(scheme_name, ["time_leap_frog_RA_space_central_diff"])
                obj.scheme_name = scheme_name;
                Cr = obj.velocity * obj.delta_t ./ obj.delta_x;

                f_llast = obj.init_func(x_grid{1});
                t_llast = obj.t_start;
                t_query_next_idx = find(obj.t_query > obj.t_start, 1, "first");

                % 若请求时刻不晚于初始时刻, 则无需计算, 直接返回初值
                if isempty(t_query_next_idx)
                    t_list(1) = t_llast;
                    f_list{1} = f_llast;

                    fprintf("%s\n\tscheme_name = %s,\n\tFinished.\n\tt_list = [%s]\n", datetime('now'), scheme_name, join(string(t_list), ", "));
                    return
                end

                if t_query_next_idx > 1
                    t_list(t_query_next_idx - 1) = t_llast;
                    f_list{t_query_next_idx - 1} = f_llast;
                end

                % 先用两层方案 (upwind) 积分一个时间步
                t_last = t_llast + obj.delta_t;
                if string(obj.bndry_func(0, 0)) == "periodic"
                    halo_left = f_llast(end-1);
                    halo_right = f_llast(2);
                else
                    halo_left = obj.bndry_func(x_grid{1}(1) - obj.delta_x(1), t_llast);
                    halo_right = obj.bndry_func(x_grid{1}(end) + obj.delta_x(1), t_llast);
                end
                f_llast_expanded = [halo_left; f_llast; halo_right];
                f_last = (1 - abs(Cr(1)))*f_llast_expanded(1 + length(halo_left) : end-length(halo_right)) ...
                        + abs(Cr(1))*(obj.velocity > 0)*f_llast_expanded(length(halo_left) : -1 + end - length(halo_right)) ...
                        + abs(Cr(1))*(obj.velocity < 0)*f_llast_expanded(2 + length(halo_left) : 1 + end - length(halo_right));
                
                % 添加首步扰动
                f_last = f_last + disturb_func_first_step(x_grid{1});
                
                % 用于记录进度
                REC_TIME_INT = 5; % seconds
                time_start = tic;
                t_now_last_rec = t_last;

                % 开始积分
                while t_last < obj.t_query(end) + obj.delta_t
                    t_now = t_last + obj.delta_t;

                    % 输出进度
                    
                    if toc(time_start) > REC_TIME_INT
                        time_remained = REC_TIME_INT * (obj.t_query(end) - t_now) / (t_now - t_now_last_rec);
                        fprintf("%s\n\tscheme_name = %s,\n\ttime_remained = %.2e,\n\tspace_dim = (%s), Cr = (%s),\n\tt_now = %.2e, t_query_next = %.2e, t_query_end = %.2e\n", datetime('now'), scheme_name, time_remained, join(string(size(X)), ", "), join(string(Cr), ", "), t_now, obj.t_query(t_query_next_idx), obj.t_query(end));
                        time_start = tic;
                        t_now_last_rec = t_now;
                    end

                    % add halo points, advance one time step
                    if scheme_name == "time_leap_frog_RA_space_central_diff"
                        if string(obj.bndry_func(0, 0)) == "periodic"
                            halo_left = f_last(end-1);
                            halo_right = f_last(2);
                        else
                            halo_left = obj.bndry_func(x_grid{1}(1) - obj.delta_x(1), t_last);
                            halo_right = obj.bndry_func(x_grid{1}(end) + obj.delta_x(1), t_last);
                        end
                        f_last_expanded = [halo_left; f_last; halo_right];
                        f_now = f_llast  ...
                                + Cr(1)*( + f_last_expanded(length(halo_left) : -1 + end-length(halo_right)) ...
                                          - f_last_expanded(2 + length(halo_left) : 1 + end-length(halo_right)) ...
                                        );
                        
                        % RA time filter
                        d_f = ra_opts.nu / 2 * (f_now - 2 * f_last + f_llast);
                        f_last = f_last + ra_opts.alpha * d_f;
                        f_now = f_now + (ra_opts.alpha - 1) * d_f;
                    end

                    % store solution at queried time.
                    if ~((obj.t_query(t_query_next_idx) - t_llast)*(obj.t_query(t_query_next_idx) - t_last) > 0)
                        if abs(obj.t_query(t_query_next_idx) - t_last) < abs(obj.t_query(t_query_next_idx) - t_llast)
                            t_list(t_query_next_idx) = t_last;
                            f_list{t_query_next_idx} = f_last;
                        else
                            t_list(t_query_next_idx) = t_llast;
                            f_list{t_query_next_idx} = f_llast;
                        end
                        t_query_next_idx = t_query_next_idx + 1;
                    end

                    % prepare for next time step
                    t_last = t_now;
                    f_llast = f_last; 
                    f_last = f_now;
                end % end of `while t_now < obj.t_query(end)`

                fprintf("%s\n\tscheme_name = %s,\n\tFinished.\n\tt_list = [%s]\n", datetime('now'), scheme_name, join(string(t_list), ", "));
                return
            end % end of `if ismember(scheme_name, ["time_leap_frog_RA_space_central_diff"])`

            warning("Invalid scheme_name. Please choose one of the following options: %s", join(AdvecProb2D.SCHEME_NAME, ", "));
        end % end of `function [f_list, x_grid, t_list] = solve_1D(obj, scheme_name)`

        function [f_list, x_grid, t_list] = solve_2D(obj, scheme_name)
            arguments (Input)
                obj
                scheme_name = "time_forward_Euler_space_upwind";
            end

            f_list = cell(size(obj.t_query));
            t_list = nan(size(obj.t_query));
            x_grid = cellfun(@(range_1D, delta_1D) (range_1D(1): delta_1D: range_1D(2)).', obj.x_range, num2cell(obj.delta_x), UniformOutput=false);
            if ismember(scheme_name, ["time_forward_Euler_space_upwind", "time_forward_Euler_space_CTU"])
                obj.scheme_name = scheme_name;

                Cr = abs(obj.velocity) * obj.delta_t ./ obj.delta_x;
                [X, Y] = ndgrid(x_grid{:});
                f_last = obj.init_func(X, Y);
                t_last = obj.t_start;
                t_now = t_last;
                t_query_next_idx = find(obj.t_query > obj.t_start, 1, "first");

                % 若请求时刻不晚于初始时刻, 则无需计算, 直接返回初值
                if isempty(t_query_next_idx)
                    t_list(1) = t_last;
                    f_list{1} = f_last;

                    fprintf("%s\n\tscheme_name = %s,\n\tFinished.\n\tt_list = [%s]\n", datetime('now'), scheme_name, join(string(t_list), ", "));
                    return
                end
                
                % 用于记录进度
                REC_TIME_INT = 5; % seconds
                time_start = tic;
                t_now_last_rec = t_last;

                while t_now < obj.t_query(end)
                    t_now = t_last + obj.delta_t;

                    % 输出进度
                    
                    if toc(time_start) > REC_TIME_INT
                        time_remained = REC_TIME_INT * (obj.t_query(end) - t_now) / (t_now - t_now_last_rec);
                        fprintf("%s\n\tscheme_name = %s,\n\ttime_remained = %.2e,\n\tspace_dim = (%s), Cr = (%s),\n\tt_now = %.2e, t_query_next = %.2e, t_query_end = %.2e\n", datetime('now'), scheme_name, time_remained, join(string(size(X)), ", "), join(string(Cr), ", "), t_now, obj.t_query(t_query_next_idx), obj.t_query(end));
                        time_start = tic;
                        t_now_last_rec = t_now;
                    end
                    
                    % add halo points, advance one time step
                    if ismember(scheme_name, ["time_forward_Euler_space_upwind", "time_forward_Euler_space_CTU"])
                        if obj.bndry_func(0, 0, 0) == "periodic"
                            halo_left = [f_last(end-1, end-1); f_last(:, end-1); f_last(2, end-1)];
                            halo_right = [f_last(end-1, 2); f_last(:, 2); f_last(2, 2)];
                            halo_upper = f_last(end-1, :);
                            halo_lower = f_last(2, :);
                        else
                            halo_left = obj.bndry_func([x_grid{1}(1) - obj.delta_x(1); x_grid{1}; x_grid{1}(end) + obj.delta_x(1)], ones(2 + length(x_grid{1}), 1)*(x_grid{2}(1) - obj.delta_x(2)), t_last);
                            halo_right = obj.bndry_func([x_grid{1}(1) - obj.delta_x(1); x_grid{1}; x_grid{1}(end) + obj.delta_x(1)], ones(2 + length(x_grid{1}), 1)*(x_grid{2}(end) + obj.delta_x(2)), t_last);
                            halo_upper = obj.bndry_func(ones(1, length(x_grid{2}))*(x_grid{1}(1) - obj.delta_x(1)), x_grid{2}.', t_last);
                            halo_lower = obj.bndry_func(ones(1, length(x_grid{2}))*(x_grid{1}(end) + obj.delta_x(1)), x_grid{2}.', t_last);
                        end
                        f_last_expanded = [halo_left, [halo_upper; f_last; halo_lower], halo_right];
                        f_now = (1 - Cr(1) - Cr(2)) * f_last_expanded(2:end-1, 2:end-1) ...
                                + Cr(1)*(obj.velocity(1) > 0)*f_last_expanded(1:end-2, 2:end-1) ...
                                + Cr(1)*(obj.velocity(1) < 0)*f_last_expanded(3:end, 2:end-1) ...
                                + Cr(2)*(obj.velocity(2) > 0)*f_last_expanded(2:end-1, 1:end-2) ...
                                + Cr(2)*(obj.velocity(2) < 0)*f_last_expanded(2:end-1, 3:end);
                        if scheme_name == "time_forward_Euler_space_CTU"
                            f_now = f_now + prod(Cr) * f_last_expanded(2:end-1, 2:end-1) ...
                                    - prod(Cr) * (obj.velocity(1) > 0)*f_last_expanded(1:end-2, 2:end-1) ...
                                    - prod(Cr) * (obj.velocity(1) < 0)*f_last_expanded(3:end, 2:end-1) ...
                                    - prod(Cr) * (obj.velocity(2) > 0)*f_last_expanded(2:end-1, 1:end-2) ...
                                    - prod(Cr) * (obj.velocity(2) < 0)*f_last_expanded(2:end-1, 3:end) ...
                                    + prod(Cr) * (obj.velocity(1) > 0 & obj.velocity(2) > 0)*f_last_expanded(1:end-2, 1:end-2) ...
                                    + prod(Cr) * (obj.velocity(1) > 0 & obj.velocity(2) < 0)*f_last_expanded(1:end-2, 3:end) ...
                                    + prod(Cr) * (obj.velocity(1) < 0 & obj.velocity(2) > 0)*f_last_expanded(3:end, 1:end-2) ...
                                    + prod(Cr) * (obj.velocity(1) < 0 & obj.velocity(2) < 0)*f_last_expanded(3:end, 3:end);
                        end
                    end

                    % store solution at queried time.
                    if ~((obj.t_query(t_query_next_idx) - t_last)*(obj.t_query(t_query_next_idx) - t_now) > 0)
                        if abs(obj.t_query(t_query_next_idx) - t_now) < abs(obj.t_query(t_query_next_idx) - t_last)
                            t_list(t_query_next_idx) = t_now;
                            f_list{t_query_next_idx} = f_now;
                        else
                            t_list(t_query_next_idx) = t_last;
                            f_list{t_query_next_idx} = f_last;
                        end
                        t_query_next_idx = t_query_next_idx + 1;
                    end

                    % prepare for next time step
                    t_last = t_now;
                    f_last = f_now;
                end % end of `while t_now < obj.t_query(end)`

                fprintf("%s\n\tscheme_name = %s,\n\tFinished.\n\tt_list = [%s]\n", datetime('now'), scheme_name, join(string(t_list), ", "));
                return
            end

            warning("Invalid scheme_name. Please choose one of the following options: %s", join(AdvecProb2D.SCHEME_NAME, ", "));
        end % end of `function [f_list, x_grid, t_list] = solve_2D(obj, scheme_name)`
    end % end of `methods (Access=private)`
end % end of `classdef AdvecProb2D < handle`

  1. Dale R. Durran (2012). Numerical Methods for Fluid Dynamics: With Applications to Geophysics [Book]. Springer New York, NY. https://doi.org/10.1007/978-1-4419-6412-0 ↩︎ ↩︎

  2. Robert, A. J. (1966). The integration of a low order spectral form of the primitive meteorological equations. J. Meteor. Soc. Japan, 44, 237–245. https://doi.org/10.2151/jmsj1965.44.5_237 ↩︎

  3. Asselin, R. (1972). Frequency filter for time integrations. Monthly Weather Review, 100, 487–490. https://doi.org/10.1175/1520-0493(1972)100<0487:FFFTI>2.3.CO;2 ↩︎

  4. Joel H. Ferziger, Milovan Perić & Robert L. Street (2020). Computational Methods for Fluid Dynamics [Book]. Springer Cham. https://doi.org/10.1007/978-3-319-99693-6 ↩︎ ↩︎

  5. Williams, P. D. (2009). A Proposed Modification to the Robert–Asselin Time Filter. Monthly Weather Review, 137(8), 2538–2546. https://doi.org/10.1175/2009MWR2724.1 ↩︎